cd /Users/Laura/Documents/Research/PDHosp/Replication
log using WMKM_REStat_state_wildci.log, replace

/************************************************************************************
PROGRAM: WMKM_REStat_state_wildci.do
PURPOSE: Calculates studentized Wild Cluster Bootstrap CIs (a la MacKinnon 2015) 
		 for state-level analysis in Wherry, Miller, Meyer, and Kaestner
DATE: 	 February 7, 2017
CONTACT: LWherry@mednet.ucla.edu
*************************************************************************************/
set seed 1028455859
global bootreps 999

***A) Studentized Wild cluster boostrap CI (does not impose null)

***weighted
cap prog drop wildst
prog def wildst

foreach y in lognp logchronic lognchronic {
	use src/`1', clear
	*ols with clustered std errors
	drop if `y'==.
	regress `y' postdiscont `2' [aweight=population], vce(cluster fips)
	global maint=_b[postdiscont]/_se[postdiscont]
	global mainb=_b[postdiscont]
	global mainse=_se[postdiscont]
	predict epshat, resid
	predict yhat, xb
	sort state
	save data/temp, replace

	*figure out number of states
	qui by state: keep if _n==1
	qui summ
	global numstates=r(N)

	postfile bskeep t_wild using $output/bs_results, replace
	*wild bootstraps
	forval b=1/$bootreps{
		use data/temp, clear
	
		*with 50% prob construct dummy that + or - Radamaker error
		qui by state: gen temp=uniform()
		qui by state: gen pos=(temp[1]<.5)
		gen wildresid=epshat*(2*pos-1)
	
		*construct new y
		gen wildy=yhat+wildresid
	
		*regress on all vars
		qui reg wildy postdiscont `2' [aweight=population], vce(cluster fips)
		local bst_wild=(_b[postdiscont]-$mainb)/_se[postdiscont]
		post bskeep (`bst_wild')
		}
	postclose bskeep

	*calculate CIs using wild t centered at original beta
	use output/bs_results, clear
	_pctile t_wild, p(2.5,97.5)
	local r1=r(r1)
	local r2=r(r2)
	local ub=$mainb-($mainse*r(r1))
	local lb=$mainb-($mainse*r(r2))
	di "SAMPLE `1'"
	di "VARIABLE `y'"
	di "Wild cluster studentized CI bootstrap for postdiscont"
	di "2.5 and 97.5 percentiles of t* distribution: (`r1',`r2')"
	di "95 percent CI is (`lb',`ub')"
	gen positive=$maint>0
	gen pos=t_wild>$maint
	gen neg=t_wild<$maint
	gen reject=positive*pos + (1-positive)*neg
	qui: su reject
	local sumreject=r(sum)
	local p_value_wild=2*`sumreject'/$bootreps
	di "p-value for postdiscont from wild bootstrap, not imposing null = `p_value_wild'"
}
end

*run for each sample and outcome under flexible and restricted models
global xvars post post#c.x_c post#c.x_c2 x_c x_c2 ///
i.fips i.bmonth
global xvars2 post post#c.x_c post#c.x_c2 x_c x_c2 ///
i.fips i.fips#c.x_c i.fips#c.x_c2 i.fips#post#c.x_c i.fips#post#c.x_c2 ///
i.bmonth

wildst bystate20091 "$xvars"
wildst bystate20091 "$xvars2"
wildst blackbystate2 "$xvars"
wildst blackbystate2 "$xvars2"
wildst nonblackbystate3 "$xvars"
wildst nonblackbystate3 "$xvars2"

wildst erbystate20091 "$xvars"
wildst erbystate20091 "$xvars2"
wildst blackerbystate2 "$xvars"
wildst blackerbystate2 "$xvars2"
wildst notblackerbystate20093 "$xvars"
wildst notblackerbystate20093 "$xvars2"

wildst bystate19991 "$xvars"
wildst bystate19991 "$xvars2"
wildst blackbystate19992 "$xvars"
wildst blackbystate19992 "$xvars2"
wildst nonblackbystate19993 "$xvars" 
wildst nonblackbystate19993 "$xvars2" 

***unweighted
cap prog drop wildst
prog def wildst

foreach y in lognp logchronic lognchronic {
	use src/`1', clear
	*ols with clustered std errors
	drop if `y'==.
	regress `y' postdiscont `2', vce(cluster fips)
	global maint=_b[postdiscont]/_se[postdiscont]
	global mainb=_b[postdiscont]
	global mainse=_se[postdiscont]
	predict epshat, resid
	predict yhat, xb
	sort state
	save data/temp, replace

	*figure out number of states
	qui by state: keep if _n==1
	qui summ
	global numstates=r(N)

	postfile bskeep t_wild using $output/bs_results, replace
	*wild bootstraps
	forval b=1/$bootreps{
		use data/temp, clear
	
		*with 50% prob construct dummy that + or - Radamaker error
		qui by state: gen temp=uniform()
		qui by state: gen pos=(temp[1]<.5)
		gen wildresid=epshat*(2*pos-1)
	
		*construct new y
		gen wildy=yhat+wildresid
	
		*regress on all vars
		qui reg wildy postdiscont `2', vce(cluster fips)
		local bst_wild=(_b[postdiscont]-$mainb)/_se[postdiscont]
		post bskeep (`bst_wild')
		}
	postclose bskeep

	*calculate CIs using wild t centered at original beta
	use output/bs_results, clear
	_pctile t_wild, p(2.5,97.5)
	local r1=r(r1)
	local r2=r(r2)
	local ub=$mainb-($mainse*r(r1))
	local lb=$mainb-($mainse*r(r2))
	di "SAMPLE `1'"
	di "VARIABLE `y'"
	di "Wild cluster studentized CI bootstrap for postdiscont"
	di "2.5 and 97.5 percentiles of t* distribution: (`r1',`r2')"
	di "95 percent CI is (`lb',`ub')"
	gen positive=$maint>0
	gen pos=t_wild>$maint
	gen neg=t_wild<$maint
	gen reject=positive*pos + (1-positive)*neg
	qui: su reject
	local sumreject=r(sum)
	local p_value_wild=2*`sumreject'/$bootreps
	di "p-value for postdiscont from wild bootstrap, not imposing null = `p_value_wild'"
}
end

wildst bystate20091 "$xvars"
wildst bystate20091 "$xvars2"
wildst blackbystate2 "$xvars"
wildst blackbystate2 "$xvars2"
wildst nonblackbystate3 "$xvars"
wildst nonblackbystate3 "$xvars2"

wildst erbystate20091 "$xvars"
wildst erbystate20091 "$xvars2"
wildst blackerbystate2 "$xvars"
wildst blackerbystate2 "$xvars2"
wildst notblackerbystate20093 "$xvars"
wildst notblackerbystate20093 "$xvars2"

wildst bystate19991 "$xvars"
wildst bystate19991 "$xvars2"
wildst blackbystate19992 "$xvars"
wildst blackbystate19992 "$xvars2"
wildst nonblackbystate19993 "$xvars" 
wildst nonblackbystate19993 "$xvars2" 


